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Abstract 

In this paper the normal collision of spherical particles is investigated. The particle interac¬ 
tion is modelled in a macroscopic way using the Hertzian contact force with additional linear 
damping. The goal of the work is to develop an efficient approximate solution of sufficient 
accuracy for this problem which can be used in soft-sphere collision models for Discrete 
Element Methods and for particle transport in viscous fluids. First, by the choice of appro¬ 
priate units, the number of governing parameters of the collision process is reduced to one, 
which is a simple combination of known material parameters as well as initial conditions. 
It provides a dimensionless parameter that characterizes all such collisions up to dynamic 
similitude. Next, a rigorous calculation of the collision time and restitution coefficient from 
the governing equations, in the form of a series expansion in this parameter is provided. 
Such a calculation based on first principles is particularly interesting from a theoretical 
perspective. Since the governing equations present some technical difficulties, the methods 
employed are also of interest from the point of view of the analytical technique. Using further 
approximations, compact expressions for the restitution coefficient and the collision time are 
then provided. These are used to implement an approximate algebraic rule for computing 
the desired stiffness and damping in the framework of the adaptive collision model (Kernpe 
& Frohlich, Journal of Fluid Mechanics, 709: 445-489, 2012). Numerical tests with binary 
as well as multiple particle collisions are reported to illustrate the accuracy of the proposed 
method and its superiority in terms of numerical efficiency. 

Keywords: Particle-laden flow, Discrete Element Method, collision modelling, Hertzian 
contact 


1. Introduction and motivation 

Particle laden flows and their numerical simulation are of major interest in a wide range 
of engineering applications as well as in fundamental research. A frequently used approach 
for the simulation of dynamic granular materials is the Discrete Element Method (DEM), 
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cf. Cundall and Strack (1979); Hoomans et al. (1996); [Poschel and Schwager| (2005). The 


linear and angular momentum balance of the particles is solved to obtain their translational 
and rotational velocity. The hydrodynamic interaction between particles is often neglected 
or fluid forces are accounted for by simple empirical correlations and the particle interac¬ 
tion modelled using macroscopic collision models of various types. The accurate numerical 
modelling of the collision process, however, is crucial for the quality of the simulation in a 
vast regime of parameters. 

Several numerical models for the collision process between particles and for the collision 


et ah, 

2007, 

2008, 

2009) 


2007, 2008, 2009). These models can be divided into two groups: hard-sphere models 


and soft-sphere models. The hard sphere approach does not aim to resolve the details of the 
collision process. Instead, large time steps are considered and the collisions are treated as 
quasi-instantaneous. The post-collisional velocities are calculated from momentum conser¬ 
vation between the states before and after surface contact. The reader is referred to ICrowel 
(2006) for details. Soft-sphere models, on the other hand, usually require an excessively small 
time step if physically realistic material parameters are matched. In the soft-sphere approach 
the motion of the particles is calculated by numerically integrating the equations of motion 
of the particles during the collision process accounting for the contact forces acting on them. 
Typical for all soft-sphere models is that very small time steps must be used to ensure that 
for reasons of stability and accuracy the step size in time is substantially smaller than the 
duration of contact. The soft-sphere contact forces are usually based on linear and non-linear 
spring damper models, reviews of which may be found in Kruggel-Emden et al. (2007, 2008 


2009). A commonly used model for time-resolved particle interactions in numerical simula¬ 


tions is a Hertzian contact force in combination with a linear damping (Lee and Herrmann 


1993 Kruggel-Emden et al., 2007). However, in contrast to linear spring models, no closed 


solution exists for this equation. One is hence forced to integrate numerically with a very 
small time step. This issue even leads some researchers to prefer linear spring models which 
are much easier to evaluate (Kruggel-Emden et ah, 2009). 

The discussion whether the linear or the non-linear approach is to be preferred seems 
unsettled in the community so far (Kruggel-Emden et al., 2009), and the present paper does 
not aim to compare these or to advocate one or the other. Instead, the mathematical prop¬ 
erties of the equation of damped Hertzian contact are discussed and an efficient engineering 
approximation is proposed so as to reduce the cost of this model. This can enhance the 
efficiency of DEMs employing the physically more realistic non-linear approach. 

Due to the time-step reduction required by many soft-sphere models, DEM practition¬ 
ers often modify the material parameters of the collision model to alleviate this problem 


(Hoomans et al. 1996) and allow the use of larger time steps. The idea is to make the 


collision process softer, hence longer in time, at the price of increased numerical overlap of 
particles during collisions. Obviously, the collision process cannot be lengthened arbitrarily. 
In particular, if the collision time becomes large compared to the free propagation kinematic 


time scale, the dissipation of kinetic energy may become anomalously small (Luding et al. 


1994). Thus, for certain dense multiphase flows, the permissible lengthening of the collision 


time may only be small. The choice of the stretching factor depends on the regime consid- 
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ered as well as the targeted accuracy and is within the responsibility of the practitioner. 
Especially for large simulations, even a small increase in step size can lead to significant 
savings in terms of computational resources and time. The results of the present work would 
be of use in the implementation of such schemes. 

Except for the linear models, the adjustment of parameters in the soft-sphere model is 
usually performed by trail and error. This is time consuming and prone to a sub-optimal 
choice of values. A current trend in the modelling of particulate flows is that DEMs are en¬ 
hanced by representations of the viscous effects of the continuous phase around the particles, 
cf. Uhlmann (2005); Kempe and Frohlich (2012b). In this framework, hard-sphere models 
are inapplicable as they cannot properly account for the coupling to the surrounding viscous 
fluid, thus introducing substantial numerical errors (Kempe and Frohlich, 2012a). Here, soft- 
sphere models are required, with the drawbacks discussed above. As a remedy, a systematic 
strategy was recently proposed to determine the parameters for a softened contact model 
(Kempe and Frohlich 2012a). It is based on imposing the duration of the contact between 
the particles during collisions according to some external constraint, such as a pre-selected 
time step. The coefficients in the model are then determined so as to maintain the exact 
restitution coefficient. This guarantees maximal physical realism under given constraints 
imposed by computational resources. The approach, termed Adaptive Collision Time Model 
(ACTM), was implemented and tested with particles in viscous fluids for single collisions 


(Kempe and Frohlich 2012a) as well as for multiple simultaneous collisions (Kempe et al. 


2014). 


Beyond that, the approach is very interesting for pure DEM without viscous fluid, as it 
provides an automated systematic approach to regularizing the collision process. Another 
substantial advantage of this approach is that the original physical values of the coefficients 
can be introduced as bounds so that the original model is obtained again in a regular limit 
when the collision time is sufficiently reduced. This provides optimal commodity for the 
user. 

It should be noted, that the full collision modelling procedure for particles in viscous 
flow proposed by Kempe and Frohlich (2012b), termed Adaptive Collision Model (ACM), 
consists not only of the ACTM but also a lubrication model and an Adaptive Tangential 
Force Model (ATFM). The latter accounts for the effects of tangential forces during surface 
contact. This issue is not addressed in the present paper, which focuses on the treatment of 
normal collisions in the ACTM paradigm only. 

In a simulation with many particles, each collision takes place with different velocities 
of the collision partners. Hence, when imposing the duration of contact and the restitution 
coefficient, one is forced to select the model coefficients for stiffness and damping for each 
collision individually. If no closed solution is available, this requires an iterative procedure as 
indeed used so far (Kempe and Frohlich 2012a). In the present paper, this is now improved 
by providing a direct solution to this problem, based on a systematically controlled approx¬ 
imation. The increased efficiency is demonstrated by suitable test cases and comparison to 
the original method. 

The paper is structured as follows. First, an exact formal solution of the equation of 
motion for a normal linarly damped Hertzian collision is derived using nonlinear transfor- 
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mations and a parametric series expansion. Then, a rigorous calculation of the collision time 
and restitution coefficient is carried out. As a next step, compact analytical approximations 
are developed. These formulae allow the direct computation of the physically relevant param¬ 
eters, i.e. collision time and restitution coefficient, from the intrinsic material parameters. 
Afterwards, the inverse problem is addressed. The artificial lengthening of the collision time, 
while preserving the restitution coefficient, requires the computation of the appropriate stiff¬ 
ness and damping. Finally, numerical tests demonstrate the accuracy and efficiency of the 
proposed algorithm, including test runs in typical engineering settings. 


2. Basic equation of the collision process 

The situation of a normal particle-wall collision is illustrated in Figure [lj Particle deforma¬ 
tion during contact is represented here by the overlap of the undeformed particle with the 
collision partner. The equation of motion governing the surface penetration ( = ((t) during 


the contact phase considered here is (Kenrpe and Frohlich 2012a) 


rripC = ~d( - k( 3/2 


( 1 ) 


with m p the mass of the particle, d the damping coefficient and k a stiffness parameter, 
the last two being material properties. The overdot represents differentiation with respect 
to time t. The second term on the right-hand side is the nonlinear restoring force originally 
derived by He rtz] (1882). The first term on the right hand side of ([Tj) corresponds to the 
damping, which is assumed to be linear. For d = 0 the behaviour is termed ideally elastic. 
The initial conditions at the beginning of the collision are ({t = 0) = 0, ((t = 0) = u m . 

Equation ([TJ) is more conveniently expressed in dimensionless form by defining new vari¬ 
ables r, z with t — rt* and ( = zu- m t *, which is tantamount to fixing the characteristic unit 
of velocity for the system as u in and choosing A as the (at present arbitrary) unit of time. 
The first and second derivatives of £ can then be expressed as 


C(t) = u in z{t ) 


and 




at) = -f z(t) 

v * 


( 2 ) 


( 3 ) 


Here and in all further cases of occurrence, the overdot denotes, when applied to a dimen¬ 
sionless quantity, differentiation with respect to the dimensionless time r, unless stated oth¬ 
erwise. Dividing by the characteristic unit of force m p u in /t* } one obtains the dimensionless 
equation 

t d kv 1 / 2 

Z + ~^~Z + = a ( 4 ) 


U\ n m r 


, 5/2 

m p t * 


Choosing the unit of time i* as 


t * 


rnt 


k 2 Uw 


( 5 ) 
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yields the following equation for z(r): 


z + 2A z + z 3/2 = 0 


with initial conditions z(t = 0) = 0, z(t = 0) = 1 and the parameter 


1 

2 rri p 


IA 

2 m p 



( 6 ) 

( 7 ) 


Following Kempe and Frohlich (2012a), the collision time T c , i.e. r c in the non-dimensional 
setting, is obtained as the strictly positive root of £, i.e. the point of time when the surface 
penetration of the particle-wall system returns to zero. The restitution coefficient is then 
conveniently defined as 


edry — 


Uont C{t — T c ) 


= Z[T = T c 


U\ 


U\ 


( 8 ) 


The latter expression makes clear that the only way the restitution coefficient may depend 
on material and other input parameters is as a function of the parameter A introduced 
above. Physically, one may interpret this as follows: All linearly damped normal Hertzian 
particle-wall collisions are, up to a characteristic parameter, similar. This characteristic 
parameter can be the restitution coefficient ea ry , which is easy to determine experimentally 
but cannot be calculated trivially from known material properties and initial conditions. 
The other option, introduced here, is to choose A as the characteristic parameter, which 
cannot be measured directly but is readily calculated. At the heart of the present work 
is the convenient analytical conversion between the two. The case A = 0 corresponds to 


the undamped case considered by Hertz (1882), which is simply a conservative system, i.e 


e^y = 1, and is integrated readily. Furthermore, Hertz (1882) gave the collision time as 


T c , o — 


VEE0 

r (S) 


25 m 2 

_P_ 

16 k 2 U\r 


( 9 ) 


with T denoting the Gamma function. Evaluation of all factors to four significant figures 
yields r c 0 = 3.218 in dimensionless form. 


3. Exact formal solution of the equation of motion 

3.1. Introduction 

The purpose of this section is to calculate rigorously from first principles, i.e. starting from 
the equation of motion (|6| and using mathematically justifiable techniques, the physical 
parameters pertaining to the characterization of collisions in the chosen model. Such results 
are not only useful for the purpose of theoretical investigation, but are also instructive 
from the perspective of the analytical technique. Since these results are, to the best of 
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our knowledge, not available in the existing literature, their derivation is presented in the 
following. 

Due to the fractional power of z in (|6]) and the initial condition z(r — 0) = 0, a naive 
power series solution for z in r is not possible. Furthermore, it is not desirable at all to have a 
series representation in t, since it is not a variable that is naturally small. On the other hand, 
especially in practical applications, very low values of e dry are unlikely to be of interest. By 
consequence, A is a parameter that may be assumed small for practical purposes, indeed, 
numerical examples show that if a lower bound of 0.4 is imposed on ed ry , A ^ 0.2 may be 
safely assumed. In addition, as remarked earlier, the case A = 0 has already been studied 
in detail and is well-understood (Hertz, 1882). Thus, in the following, the solution of the 


equation of motion is most conveniently expressed using series expansions in A. 


3.2. Solving the equation of motion 

It will be found that the relevant mathematical operations are substantially facilitated when 
the problem is transformed to phase space (z,z), as visualised in Figure [2j With this ap¬ 
proach, the dimensionless velocity can be interpreted as a function of the penetration, i.e. 
z = z(z). Using the relation z = di/dr = dz/dr dz/dz = z dz/dz the equation of motion 
expressed in these new variables reads 

+ 2Ai + z z/2 = 0. (10) 


In passing, observe that if the damping term were absent, this would reduce to an exact 
differential, with the solution given by the level curve of a potential function which is readily 
interpreted as the conserved total energy of the system. In fact, the energy method was the 
procedure followed by Hertz (1882) when deriving the solution for the undamped case. The 
analysis here, however, needs to be more involved. It is now convenient to define 


z =: 



in 

out 


( 11 ) 


where ‘in’ refers to the part of the collision during which the particle is compressed and the 
motion is directed into the wall, while ‘out’ describes the outward motion of the particle 
(Figure^. By the chain rule, one has zdz = ldu±, where the choice of + and — is determined 
by 0- Finally, it is helpful to define y := y/z, whence dz = 2ydy. Inserting these new 
variables into (10) yields the initial value problems of the two phases of motion: 


-r^ + 8A yy/vf + 4r/ 4 = 0 
_ 

zyj- - 8A y^/vZ + Ay 4 = 0 


v+(y = 0) = 1 
v-(y = 0) = e 2 dry . 


( 12 ) 

(13) 
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For A assumed small, the unknowns are now expressed in form of a series in powers of A: 


V+ = U+,mA r 

m =0 
oo 

Y 


V _ = 


^+,m 


u_ ™ = 


1 <9 m u + 


m\ dX m 
1 d m v _ 


m=0 


m! <9A m 


A=0 


A=0 


(14) 

(15) 


Differentiating m times with respect to A and evaluating at A = 0 successively yields the 
equations determining v±_ m for m ^ 0. For m — 0, this gives 


4„,5 


and for m ^ 1, 


v±, 0 = 1 - g y 


f) v , P) m 

m\ —— = T% 


A=0 


<9y <9A r 

The right-hand side of (JTt]) may be further simplified using the Leibniz formula, 


Qr. 


dX> 




£ 


m\ d l X 
dX< 




A=0 


dX r 


(16) 

(17) 

(18) 


A=0 


A =0 i=0 

Only the term i = 1 survives, so that the higher order corrections are given by 

' y d m “ 1 DF7 


^+,m 


u_ ™ = 


m 


1) 


if 


1 

r) m P 2 
° e dry 

m\ 

3X m 


+ 


5 A ”"- 1 
8 


y'dy'; 


A=0 


A=0 


[m 


1) 


(I 




<9A m_1 


y'dy'. 


(19) 

( 20 ) 


A=0 


Since the (m — l)-st derivative of v± at A = 0 can only involve u± i0 , • • •, equations 

(19) and (|20| yield well-defined recursions for v±^ m . Furthermore, it is readily seen by math¬ 


ematical induction, that v± tTn have a convergent Taylor series expansion around y — 0, and 
that the convergence radius is (5/4) 1 ^ 5 . From conservation of energy, however, it can be 
deduced that z ^ (5/4) 2//5 , i.e. v± tTn can be represented by a power series on the whole 
domain of interest. A physical interpretation of these terms is discussed in |Appendix A 


Note also that derivatives of e dry appear explicitly in the expressions for u_ jm . This does not 
imply that e dry needs to be known a priori. Rather, it will be seen later that e dry is uniquely 
determined by certain conditions that need to be fulfilled in order for the trajectory to be 
continuous in time and space. 

To find a connection to the time domain, one may exploit the fact that z dr = dz to find 


T. 


T = T (y) = 


+ 


(y) := f 

Jo 


2d y'y' 


m 


T-(y) ■= T+(y) + f 

Jo 
7 


2dy'y' 


( 21 ) 
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with y 2 = z the maximum surface penetration. This solves the problem in the classical 
sense. In order to obtain an explicit relation of the form z = z(t), one would have to 
take the functional inverse and square the resultant expression. However, since it does not 
contribute to the calculation of the relevant physical properties z, e^y and r c , that part of 
the problem, though utterly non-trivial, will not be addressed here. 


3.3. Calculation of the physical parameters 

With the above solution (19)—(21) of the equation of motion at hand, it is now possible to 
calculate the significant physical quantities. Although the two primary parameters of interest 
are the restitution coefficient ejry and collision time r c , the calculation of the maximum 
penetration £ is required as an intermediate step and is therefore addressed first. From the 
definition of z, it follows that z = 0, which is equivalent to the more tractable condition 


v+{y = V) = 0 , (22) 

with y 2 = z the maximum penetration. This equation uniquely determines z in terms of the 
smallest positive real root of v + . Since the solution trajectory is required to be continuous, 
the velocity on both the inward and outward trajectory must go to zero at the same value 
of y, yielding the condition 

v-(y = y) = 0. (23) 

This uniquely determines e dr y Finally, the collision time can be computed directly as 

T c = r-(y) (24) 


Obviously, each of the above quantities depend on A. For A = 0, the classical case of the 
undamped Hertzian restoring force is obtained, for which the solutions are (Hertz, 1882) 


2/o := ?/(A = 0) = (f) 1/5 , (25) 

e dr y(A = 0) = 1 , (26) 

t c ,o := r c (A = 0) = 3.218 . (27) 


For A small, therefore, it is natural to develop the solutions in terms of a power series in 
A around A = 0, whereby the respective governing equations are successively differentiated 
at A = 0 using implicit differentiation, with the chain rule used to evaluate the Taylor 
coefficients. 


3-4- Explicit computation of first-order corrections 

Following the procedure specified above, it should, in principle, be possible to calculate 
the correction terms to arbitrarily high order, which naturally becomes more and more 
involved with increasing order. Here, the explicit computation of the first-order corrections 
is presented for illustration. 

Maximum surface penetration To zeroth order, as established above, condition (22) yields 
the result y = y/5/4 + 0( A). This is the same solution as one would estimate using conser¬ 
vation of energy, since to zeroth order in A, the dissipation of energy, i.e. work done by the 
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damping force is negligible. Differentiating v + (y) = 0 once with respect to A and using the 
chain rule to account for the fact that y = y(A), one arrives upon evaluation at A = 0 the 
condition 


dv+,o 

dy 


d y 


. 8X 


y=vo 


+ v+,i(y = y 0 ) 

A=0 


0 , 


or 



dy 

d\ 


A=0 



| y 5 y d y 


0 


(28) 


Carrying out the integration yields 


V = 



1 Jwyr® 

2 V 25 mRm) 


A+0(A 2 ) . 


(29) 


Restitution coefficient Recall that imposing continuity on the solution yields the condition 
V-(y) = v + (y) = 0. Again, to zeroth order one finds ed ry | A=0 = 1, which is an obvious 
consequence of conservation of energy (the dissipation through damping being first order 
and higher in A). Differentiating the expression for v_(y) yields 


dv- 


,o 


dy 


y=y o 


dy 0 

dX 


+ v-p. (y = y 0 ) = 0 


A=0 


and hence 


de, 


dry 


d X 


4 dy 


A=0 


= 2S » A 


A=0 ^0 


^yo j - rv o j - 

4 Jl - fz/ 5 y dy = -8 J 1 - ±y 5 y dy . 

J 0 v J 0 v 


(30) 


Carrying out the integration yields 


e dr y — 1 - 2 \ 7 X + C>(A 2 ) 


16 Arf 1 ) 

io Viol 


(31) 


Collision time The collision time requires a computation of the integral 


Tr = 


I 


y / 1 1 
+ 


2 ydy 


(32) 


Since v± = u ± i0 + u±,iA + •••■, this gives 


Tr = 


f 

f 

Jo 


- 1/2 


—= (1 + ^A + G(A 2 )^) +—— (l + — A + CXA 2 1 

y/V+fi V U +.0 J V V ~,0 


- 1/2 


2 y dy 


+ 


/25vAFr(|) 2A 


1-h 5 Vl 6 D r (D) 


1 - # 2/ 5 


+ 0{ A 2 


2 y dy 


(33) 
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using (16), (19), (20) and (31), as well as the Taylor series representation of (1 + x) 1 X 1 . For 


A = 0, the integral is straightforward, yielding the same result as obtained by Hertz (1882) 


V? r (D 



(34) 


For the first-order correction, differentiating with respect to A using differentiation under 
the integral sign (Leibniz rule) and setting A = 0 gives 


dr c 


dX 


\=o 



+ 


y=vo 


5 /25 y^T(|) p Aydy 

V 16 ^ r (^) Jo ^1- | ys'j 


(35) 


Both terms are divergent. The first term on the right hand side has a simple pole at y — y 0 
and the integrand in the second term diverges too fast for y —y y 0 . The result can, however, 
be regularized by calculating the integral first for general y < y 0 and then taking the limit 
of the whole expression for y —* y 0 , i.e.: 



Evaluating the integral in the second term of the right-hand side yields a contribution 
~ y 2 ( 1 — |y 5 ) -1 and a second one involving a hypergeometric function. One arrives at the 
following expression: 


d r c 
~dX 


A=0 


v^r(j) 
ads) 


= lim 


1 2 y ( y/y 0 - 1) 


vl 


1 - (yfy o) 


+ 2 4v 2 ^( hbb(y/y°r 


(37) 


where 2 F \(•) denotes the (ordinary) hypergeometric function (cf. Abramowitz and Stegun 


1965, Chapter 15, pp. 555). Note that, wherever opportune, the factors have been expressed 
in terms of y 0 for ease of manipulation. The first term in the limit can be evaluated using 
de l’Hopital’s rule and is found to vanish for y —>• y 0 . The second term remains regular in 
the limit and can be simply evaluated at y = y 0 , yielding 


dr c 

d\ 


A=0 


= 2 4 v^£(i) 

c y 0 9 r/ 9 \ 

0 io 1 VioJ 


=n(f,|;|a) ■ 


(38) 


This can be simplified using Gamma functions, so that one finally has 


T r = 


2 v vr(j) s /25 , ^5idr(!)' 


F (S) 


l 6 + 


r(i) 


A + 0 (A 2 ) . 


(39) 
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3.5. Summary of analytical results 

To summarize, the upshot of the theoretical calculations presented in Section [3] are expan¬ 
sions of the relevant physical quantities in terms of the (small) parameter A. Thus, for the 
square root of the maximum penetration 

\Tz = y 0 - 0.748A + 0.578A 2 + C>(A 3 ) . (40) 

This finally yields: 

e dry = 1 - 3.576A - 5.131A 2 + C>(A 3 ) (41) 

t c = t c ,o + 1.152A + 0(A 2 ) (42) 

where the special functions and other constant numerical factors were evaluated to four 
significant figures. 


4. Approximate method and results 


While the solution (41), (42) derived in the previous section is formally correct, calculating 


the necessary quantities can be costly as higher-order terms must be included for sufficient 
accuracy. If, however, terms higher than quadratic are involved the expressions cannot be 


conveniently inverted. Indeed, the inversion of (41) and (42) cannot be carried out by al 


gebraic means if quintic or higher-order terms are present. The inversion could be achieved 
using the theta function, but this is a non-elementary special function and costly to eval¬ 
uate. In the following, a more elegant procedure is developed, based on the fact that the 
expressions of e dry (A) or r c (A) can be resolved for A. 

To this end, consider the following ansatz with A, B,C denoting free parameters whose 
values are to be determined appropriately: 

aAr c0 


e dry = exp 


T r = 


y/l-CX 

t c ,o 


y/l- AX - B A 2 


(43) 

(44) 


The reason one may expect (43), (44) to be a more convenient choice for describing the 


required quantities than the simple power series in (41) and (42) is because they mimic to a 


great extent the relationships that are valid in the case of the damped harmonic oscillator 
which is obtained if the power | is replaced by unity in (|6]). The equation of motion of the 
harmonic oscillator is readily solved in closed form and yields the behaviour loge dry oc — 5t c 
and r c oc (1 — <5 2 ) -1 / 2 , where 5 corresponds to the damping ratio. 


The value of most of the free parameters in (43) and (44) can be determined directly 
from the solution (41) and (42) found in Section |3j Differentiating (43) at A = 0 yields the 
two equations 


de 


dry 


d\ 

d 2 e 


QC7"c,0 j 


A=0 


dry 


dX 2 


~c,0 (^^”c,0 


(45) 

(46) 


A=0 
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The derivatives appearing on the left-hand side are given by (41), whence one obtains 


a = 1.111 and C = 0.744 


when evaluated to four significant figures. Likewise, differentiating (44) yields 


dr c 

dX 


= \t c ,qA , 


A=0 


so that 


A = 0.716 


( 47 ) 


(48) 


(49) 


The method is difficult to employ for determining the remaining free parameter R, because 
it requires knowledge of the second-order correction to r c , which in turn necessitates the 
summation of divergent terms (in the form of derivatives of terms like (1 — | y 5 ) -1 / 2 with 
respect to y at y = y 0 ) and improper divergent integrals. This appears, in fact, to be a char¬ 
acteristic shared by corrections to the collision time of all non-zero orders. For the first-order 
correction, it was possible to carry out the regularization of the result by hand, because the 
indefinite integral underlying the divergent improper integral could be formulated in closed 
form and the cancellation of the divergent terms could be done away with in a straightfor¬ 
ward manner. For higher order corrections, the integrals become increasingly involved; so 
much so, that even for the second-order term, the integration could not be carried out by 
hand. The calculation of B by self-consistent methods is, hence, a highly non-trivial task 
which may be addressed in a future investigation. 

For the present study, it shall suffice to obtain an approximation by a method similar 
to data-fitting. To this end, the dimensionless collision time is determined from numerically 
calculated solutions of the governing equation (J6]) for several values of A. The numerical 
integration was carried out using a Runge-Kutte method with an adaptive step size, em¬ 


ploying the routine ode45 of MATLAB™. Equation (44) predicts that (t C)0 /t c ) + AX is a 
linear function of A 2 . Linear regression of the same yields 


B = 0.830 


(50) 


The result of this section now is the approximation (43), (44) with the parameter val¬ 


ues given in (47), (49) and (50). This shall serve as a convenient alternative to the series 


expansions (41), (42) for the purposes of the subsequent sections. 


In spite of the simple nature of the approximate formulae proposed in this section, their 
accuracy is very high in the practically interesting parameter range (A ^ 0.2). For a given 


value of A, the values of ed ry and r c calculated from the closed-form expressions (43) and 


(44) were compared with the corresponding values obtained from a numerical Runge-Kutta 


solution of equation Q. It was found that the expression (44) for r c has a maximum relative 
error of 2.2 x 10 -4 . Likewise, the relative error of the approximation (43) does not exceed 
4.9 x 10 -3 for the values of A considered here. This may be seen as a numerical validation 
of the solutions derived in Section 13.31 
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5. Application to the inverse problem 


5.1. Inverse problem 

In practice, one is often interested in a sort of inverse problem to the one discussed above. The 


ACTM of Kernpe and Frohlich (2012a) imposes e dry and T c and determines the values of the 


parameters k and d such that this is obtained individually for each collision. The problem, 
then, is to chose the appropriate parameters k, d depending on the incident velocity w in , given 
the mass of the particle m p , collision time T c and restitution coefficient e dry . To the best of 
our knowledge, a non-numerical means of achieving this is not known in the literature, so 


that Kernpe and Frohlich (2012a) resorted to an iterative numerical method employing a 


Newton-type scheme briefly recalled below. 

5.2. Direct determination of damping and stiffness 


Equation (43) allows to solve for A, yielding 


A = 


« 2 t- c 2 o 


-\Cr) + yj\C 2 rj 2 + a 2 rf 0 rj 


(51) 


where r/ = (loge dry ) 2 for convenience. Only the positive root of the quadratic equation 


is physically relevant. The values to be used in (51) are r Cj0 = 3.218 (Hertz, 1882) and 
a = 1.111, C = 0.744 from Section |4} When the collision time T c is imposed in physical 
units, relation (43) can be inserted into (|44j) in order to calculate the unit of time 


U = — Vl- AX- BX 2 , 

T c,0 


(52) 


where the values A = 0.716 and B = 0.830 from Section H] are to be used. The definitions 
of t* and A, ([T]) and ([5]) respectively, then yield the physical material parameters 


, 2A m v 
d = 


k = 


t * 

m r 


y/ ^in 


(53) 

(54) 


For given material parameters of the particles, m p and ed ry , this provides a direct, non¬ 
iterative method to compute the parameters k and d in the governing equation (JTj) from 
given impact velocity u [n and desired collision time T c . 


6. Numerical results 

6.1. Validation and assessment for binary particle-particle collisions 

In this section numerical tests of the proposed subroutine for parameter adaption are pre¬ 
sented. In particular, the accuracy and efficiency of the results obtained using the developed 
method is investigated. Table [l] provides data of test runs for various values of e dry and com¬ 


pares the performance with that of the quasi-Newton iterations developed by Kernpe and 


[ 


] 
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Frohlich (2012a) to compute the desired damping and stiffness. Since the standard Newton 


iteration scheme requires the calculation of the Jacobi matrix for each iteration step, it was 
replaced by a Broyden approximation in that reference. The convergence criterion in the 
iterative scheme was based on the residual of the rebound velocity and collision time and 
was set to 10 -6 . A striking feature of the new algorithm is that the computation time is 
practically negligible for all cases. This is because all the underlying calculations involve at 
most the evaluation of elementary mathematical functions. On the other hand, the CPU time 
used by the iterative procedure shows a general trend of increasing CPU time for decreasing 
values of e dry , and always yields larger times than that required by the present method. 

Moreover, the present subroutine does not trade off accuracy for efficiency, at least not to 
any significant extent. To quantify this aspect, the equations of motion were solved with the 
adapted parameter values k appT , d appr given by the new method using a third-order Runge- 
Kutta solver. This is similar to how the algorithm would be employed in practice, as shown 
in the next section. On this basis the restitution coefficient e d p y r and the collision time T c appr 
were determined. The quality of the result can then be quantified by the relative errors 


c appr 
e dry ' 


appr 

1 C dry 



■Tcappr 

and 

appr_ 

e T c • — 

i 1 c 

^dry 

T 


(55) 


with respect to the preset target values e dry and T c . From Table [TJ it is apparent that accu¬ 
racy is indeed not compromised by the new direct method. In all cases of potential practical 
relevance (ed ry > 0.7) the relative error in the restitution coefficient is e appr ~ 10 -4 and 
does not exceed 1.3%o, which is excellent. For obvious reasons, the accuracy decreases when 
lowering e dry . Even for e dry as low as 0.4, the relative error is e app J ~ 3%, which is acceptable. 
In the time domain, the collision time that results from the parameter adjustment by the 
present method agrees up to a relative error e^? pr ~ 10 -4 with the pre-set target value. 

The computation times for the Newton iteration scheme given in Table [l] were obtained 
for a much more stringent convergence criterion (10 -6 , as stated above). Hence, one may ask 
what the computational effort is, if the required accuracy of the Newton scheme is reduced 
to that of the direct approximate scheme, simply by stopping the iterations earlier. The 
resulting timings are reported in Table [2] It turns out that in the extreme case of e dry = 0.4, 
the CPU time required when relaxing accuracy is reduced by a factor of nearly 6.3. For 
e dry = 0.95 used later in Section 16.21 the CPU time is practically not reduced at all when 


lowering the accuracy of the iterative scheme. In all cases, even with the relaxed convergence 
criterion, the required time is orders of magnitude higher than the direct scheme, because 
even the initialisation and execution of 4-5 iteration steps is much costlier than the evaluation 
of elementary functions. 

In summary, it may be concluded that the method proposed in this paper enables a 
computationally cheap and at the same time sufficiently accurate implementation of the 


concept of parameter-adapted time-stretched collision modelling developed by Kempe and 


Frohlich (2012a) 
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6.2. Application to multiple particles sedimenting on a rough surface 

The final test case reported here deals with the sedimentation of 100 randomly distributed 
particles and their impact on a layer of 195 fixed particles arranged in hexagonal packing 
(Fig. [3^) very similar to the simulations presented by Kempe et al. (2014). The main dif¬ 
ference is that the fluid is neglected in the present case corresponding to an infinite Stokes 
number of the particles. The goal of the simulations is to compare the results obtained with 
the original scheme of Kempe and Frohlich (2012a) when applied to this case with infinite 
Stokes number, to the results of the approximate inverse procedure developed in Section [5] 
The situation is very close to practical applications as it features multiple simultaneous 
collisions. At this point, it should be recalled that, in keeping with the pre-declared scope 
of the present work, the simulation presented below only considers normal collision forces, 
since only these enter the computations of the ACTM. In other words, the collisions are 
generally oblique, but with the tangential forces supposed to be negligibly small. The goal 
of the test is to demonstrate sufficiently close agreement of the results obtained using the 
direct algorithm with those of the original iterative procedure and futhermore to obtain 
timings for such a realistic situation. 

The computational domain is = [0, L] x [0,L] x [0, L] with L = 1.5. Periodic boundary 
conditions were applied in the x- and ^-direction. The gravitational acceleration is g = 9.81 
and the density of the particles is g p = 1200, with the particle diameter D = 0.1154. The 
mobile particles are placed randomly in the subdomain fA = [0.1,1.4] x [0.3,1.2] x [0.1,1.4], 
and their velocity initialized with zero. Once these particles are released from their initial 
position (Figure^), they are accelerated by gravity towards the fixed layer and then collide 
with the bed or with other mobile particles. At the same time they are subjected to a 
dissipation of kinetic energy when undergoing collisions. Two different cases are considered 
here. In Case 1, the coefficient of restitution is e^y = 0.95 and in Case 2 the value is 


Cdry = 0.7. The simulations were run with At = 5 x 10 -4 for 5000 steps, corresponding to a 
non-dimensional simulation time of t = 2.5 . This situation is depicted in Figure^ for Case 
1. The particles are coloured according to the absolute value of their velocity u p = |u p | from 
red (tip = 2) to blue {u p = 0) showing that not all the particles come to rest at the end of 
the simulation if the damping is weak. 

The collision process is elucidated by computing the various components of the total 
energy of the particles. The potential energy, the kinetic energy and the energy stored by 
deformation of the springs in the collision model are defined as 


and 


n p 

E pot ^ ^ WTp,i 9 Vp,i 

i= 1 



i =1 



(56) 

(57) 


( 58 ) 
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respectively. The total energy in the computational domain then is 

£tot £pot T -^kin "T -^spring (59) 

The fractions of energy are displayed in Figure [4}a and Figure [4jo for Case 1 and 2, 
respectively. As already mentioned above, in Case 1 the particles do not come to rest at 
the end of the simulation, reflected by their kinetic energy not being zero at the end of the 
simulation. In contrast to this, the kinetic energy of the particles at t — 2.5 is zero for Case 
2 . Obviously, no significant differences of the various fractions of the energy are observed if 
the iterative or if the approximate method is used. The slight differences that are seen in 
Case 1 rather result from the fact that extremely small differences in the collision process 
can yield different particle trajectories, in turn leading to subsequent differences, as in a 
billiard system, without changing the overall bulk behaviour of the ensemble. 

To further elucidate the efficiency of the two numerical procedures, the CPU times of the 
old and the new scheme are compiled in Table [3] Here, t tot is the overall CPU time required 
for the whole simulation including overhead, t par the time spent in the particle routines, t coe ff 
is the part of t par required for the determination of stiffness and damping, n co i the overall 
number of collisions and, finally, t co i is the average time required per collision. Obviously, 
the numerical effort is substantially reduced for all values of e dry if the new approximate 
method is used. The results presented in this section confirm the accuracy, robustness and 
efficiency of the proposed method. 


7. Concluding remarks 


In this paper, normal particle-wall and particle-particle collisions were studied by modelling 
the elastic interaction with a repulsive Hertzian contact force and an additional damping 
force linear in the velocity. 

First, the equation of motion ([l]) was converted to its dimensionless form ([6]). This re¬ 
duces the number of parameters in the equation to a single constant A depending on the 
material parameters k, d and the impact velocity w in . In particular, the physically relevant 
and practically interesting properties, the collision time T c (dimensionless: r c ) and restitution 
coefficient e dry , hence, only depend on this particular combination of physical input param¬ 
eters. This is conceptually and technically analogous to the damping ratio, or its inverse, 
the quality factor, which are familiar from the simple harmonic oscillator, the latter more 
so from LC circuits in electrical engineering but appears to be new in the present setting of 


damped Hertzian collisions. In contrast, earlier works (Kruggel-Emden et al. 2007; Kempe 


and Frohlich 2012a) typically use a three-parameter family of input variables k , d, u in . While 


the use of e dry to label cases is helpful in practice, the proposed nondimensionalisation and 
the subsequent analytical investigation establishes in a straightforward manner the one-to- 
one correspondence between the important experimental parameter e dry and the parameter 
A, which is easy to obtain from material parameters. Furthermore, it demonstrates that dif¬ 
ferent parameters with the same e dry have solutions that are identical up to a linear scaling 
of space and time. 
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Next, an exact formal solution of the governing equation using nonlinear transformations 
and a series in the parameter A was proposed, providing a rigorous calculation of e dry and r c 
from the equation of motion. Owing to the technical difficulties presented by the governing 
equation, the analysis is necessarily involved, the methods and technique employed may be 
instructive, and the results may be of use in theoretical considerations. The approach is also 
quite flexible, and a similar analysis would be applicable to more general settings, such as 
the extension to fully nonlinear models, i.e. a nonlinear spring force in combination with a 


nonlinear dissipative force as considered in some studies (Kawabara and Kono, 1987 Crowe 


2006 Crowe et al. 1998). This may be the subject of future investigation. 


Subsequently, an approximation based on the exact solution was proposed, yielding com¬ 
pact and convenient formulae for e dry and r c . Inverse formulae were derived on that basis, 
enabling a very efficient and accurate calculation of the required stiffness and damping from 
the given input parameters, i.e. the desired collision time T c and restitution coefficient e dry . 
These were then applied to an engineering context. The quasi-Newton iteration scheme in 
the original ACTM was replaced by the direct approximate solution developed here. Nu¬ 
merical tests on binary and multiple particle collisions confirm the accuracy and efficiency 
of the proposed method. 
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Appendix A. Physical interpretation of the coefficients v± jm 

It may be instructive to discuss the physical meaning of the coefficients v± >m in the series 


expansion in powers of A of v±, equations (19)—(20). First, it is readily seen that the zerotli- 
order term is 


4„,5 


v+,o = o = 1 - Yj y 


as obtained in (16). The fact that the two expressions are equal indicates reversibility. Indeed, 


for A = 0, the classical oscillation problem with Hertzian restoring force is obtained, which is 
a conservative system. Moreover, upon inspection of the functional dependence, it is evident 
that the zeroth-order term expresses conservation of mechanical energy itself. The kinetic 
energy per unit mass and in dimensionless form is given by 

T/kin 2^ 5 

depending on which part of the trajectory is considered. Thus, the auxiliary variables v± are 
proportional to the kinetic energy. The potential energy associated with the Hertz contact 
force (dimensionless and per unit mass) is 

£ pot = §= §y‘ . 

For both the inward and outward motion of the particle, the conservation of energy 

\v± + | y 5 = E kin + E pot = const. (A.l) 

is valid to zeroth-order in A. 


Now consider the more involved first-order corrections. Equations (19)—(20) yield 


v+,i = 


v -,i = 


f 


i-ly 5 y'dy' 


d e l r y 


d\ 


+ 


A=0 


8 fo ' 


(A.2) 

(A.3) 


Both integrals can be evaluated analytically to give 


„ _ 16 2 E ITT 20 „ 2 Tp (2 1.7. 4 5\ 

v +,i - -^y \f l ~- 5 y ~ ~9~y 2 ^H5> 2> 5> 5 y) 


with 2 F\ denoting the hypergeometric function (cf. Abramowitz and Stegun 1965, Chapter 
15, pp. 555). A similar expression results for U-p. The integral representations are actually 
more conducive for physical interpretation. Indeed, both expressions can be interpreted as 
work done by the dissipative force. In case of u_ a , there is an added constant term, because 
the initial condition (y = 0) = e^ ry depends on A, unlike the much more straightforward 
case v + (y = 0) = 1. For this reason, the +, or inward, branch of the trajectory shall be 
discussed in the following. Since 2 ydy = dz , 


Ki=-2 rvi-c /2d ‘'- 


(A.4) 
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On the other hand, from v + = n +j0 + 0( A), it follows that to first order in A, the work done 

by the damping force is 

z d z' = —2A ^ y/v±fi dz! + 

Jo 

Thus, v± :1 are essentially first-order corrections in the energy balance due to the work done 

by the dissipative term in the equations of motion. 
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Figure 1: Sketch of collision problem with particle position over time 




Figure 2: Sketch of the compression (in) and the rebound phase (out) during the collision 
process (cf. Figure [l| with 5 being the non-dimensional maximum surface penetration, a) 
Physical space, b) Phase space 
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Figure 3: Sedimentation of 100 particles at infinite Stokes number and impact on a layer of 
195 fixed particles in hexagonal packing with e^ry = 0.95 (Case 1). a) Initial configuration 
at t — 0, b) situation at the end of the simulation t = 2.5. The particles are coloured by the 
absolute value of their velocity from red (u p = 2) to blue (u p = 0). 




Figure 4: Fractions of energy for the sedimentation of 100 mobile particle and collision with 
a layer of fixed particles using the ACTM (Kempe and Frohlich 2012a), labelled ‘iterative’, 
and the present approximate scheme (51)- (54). (a) Case 1 with e^ y = 0.95, (b) Case 2 with 
e dry = 0.70. 
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^dry 

< 

^CPU [ S ] 

l 

— 

— 

0.95 

4 

0.0234 

0.90 

4 

0.0234 

0.80 

4 

0.0234 

0.70 

4 

0.0234 

0.60 

4 

0.0234 

0.50 

5 

0.0313 

0.40 

5 

0.0313 


Table 2: Number of iterations n[ t and computation time t' CPlJ for the quasi-Newton itera¬ 
tion scheme of Kempe and Frohlich (2012a) with the convergence criterion at par with the 
accuracy achieved by the present method (cf. Table [lj). 


Method 

e dry 

^tot [s] 

^par [s] 

^coeff [s] 

^col 

^col [s] 

Newton 

0.95 

836.46 

812.78 

221.49 

21950 

0.01009 


0.90 

839.22 

815.23 

229.71 

22362 

0.01027 


0.80 

783.91 

760.36 

186.57 

16605 

0.01134 


0.70 

734.28 

711.93 

142.86 

11687 

0.01222 

present 

0.95 

597.76 

578.45 

1.953 x 10“ 2 

21111 

9.72 x 10“ 7 


0.90 

601.75 

582.09 

3.516 x 10“ 2 

22974 

9.05 x 10“ 7 


0.80 

591.66 

572.28 

2.344 x 10“ 2 

18845 

8.32 x 10“ 7 


0.70 

591.19 

571.71 

7.813 x 10“ 3 

11829 

8.38 x 10“ 7 


Table 3: CPU time required for the sedimentation test case with various values of the resti¬ 
tution coefficient using the ACTM with Newton iterations (Kempe and Frohlich, 2012a) and 
the present direct method (51)—(54) for the computation of stiffness and damping. Nomen¬ 


clature defined in the text. 
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